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ABSTRACT 



Chaos describes a class of motions of a deterministic system whose time 
history is sensitive to initial conditions. Because of the sensitivity of initial 
conditions, the response of a dynamical system may result in instabilities. 
Hence, a study of nonlinear response of structures under the expected 
frequencies of excitation becomes important. Chaotic behavior, for example, 
may be found in the vibration response of large flexible space structures 
including trusses, booms, and radio antennas. Methods of quantifying chaos 
have been applied to flexible beams both analytically and experimentally. This 
research effort investigates the effects of sensors, strain gages and 
accelerometers, in studying chaotic motions. A long flexible beam is used to 
model the chaotic behavior, which is also mathematically modeled as Duffing's 
Equation. Time histories are recorded and analyzed using pseudo-phase space, 
Fourier spectrums, Poincare sections, Lyapunov exponents and fractal 
correlation dimensions. Comparison of the two sensors is also performed. 
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I. INTRODUCTION 



A. REASON FOR ANALYSIS 

Many of the satellites currently being launched have highly flexible 
components such as trusses and booms. Long, highly flexible booms and 
trusses are used for gravity gradient stabilization, supporting large solar arrays 
and radio antennas, and ensuring that humans and equipment are adequately 
shielded from radioactive power generators in some proposed interplanetary 
spacecraft. Deploying these trusses and booms and maneuvering the satellite 
for attitude control and station keeping result in large displacements and 
rotations of the tips. This leads to the presence of geometric nonlinearities in 
the structure. Material nonlinearities are not usually present since the booms 
and trusses are not deflected past their yield point and thus the linear stress- 
strain law applies. Because of the geometric nonlinearities, it is appropriate to 
model flexible space structures by nonlinear equations of motions instead of 
linear ones. Duffing’s Equation is one equation used to adequately describe the 
nonlinear behavior due to the elastic effects of these large structures. 

The presence of nonlinearities may result in motions which are impossible 
to predict after long periods of time, yet completely deterministic since the 
equations of motions are known. Slight changes in initial conditions produce 
radically different outputs. These types of responses have been characterized 
as "chaotic." Chaotic vibrations of spacecraft structures could result in control 
problems. 
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B. SCOPE OF THESIS 

The purpose of this theses is to demonstrate that chaotic vibrations can 
occur in highly flexible space structures using theoretical and experimental 
analyses. For the theoretical analysis, a flexible sheet of aluminum 7 ft. long, 3 
in. wide, and 0.080 in. thick, was modelled by Duffing's Equations. This process 
is described in Chapter III. For the experimental analysis, the same beam was 
instrumented with strain gauges and accelerometers. The beam was then 
shaken by an MB Dynamics Model PM500A Vibration Exciter and the time 
histories of the strains and accelerations were recorded. This process is 
discussed in Chapter IV. Data from both analyses were then analyzed using 
the computer program, "Chaos," developed by CDR Martinus Sarigul-Klijn, 
USN, [Ref. 1]. 

In addition, the purpose of this thesis is to investigate the sensitivity of 
sensors in detecting chaotic vibrations. The two most commonly used types of 
sensors are accelerometers and strain gauges. Measuring accelerations may 
be a more accurate method of recording the time histories of chaotic systems 
since this includes inertia, damping, and stiffness terms of the mathematical 
equations used to model the system. Strain gauges permit measurement of only 
the stiffness terms of the mathematical equations used to model the system. 
After a discussion of the Science of Chaos (Chapter II) the results of the 
theoretical and experimental analyses will be discussed. 
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II. THE SCIENCE OF CHAOS 



A. BACKGROUND AND PREVIOUS STUDIES 

A deterministic system whose time history has a sensitive dependence on 
initial conditions is one definition of a chaotic system [Ref. 2: p. 4]. Even though 
the equations of motions which adequately describe a system are well known, 
small changes in the initial conditions produce markedly different outputs. For 
example, although there are many computer programs which can accurately 
predict the weather tomorrow, there are no programs which can accurately 
predict the weather for a day next year. Recognizing that chaos occurs in the 
nonlinear deterministic models raises the hope of understanding the source of 
this random-like behavior and possibly doing something to change it. 

Although Henri Poincare observed chaotic behavior in the trajectories of 
some celestial bodies at the beginning of the twentieth century, the science of 
chaos is a relatively new field. It is one which essentially became possible after 
the birth of the computer which not only allowed numerical integration of 
systems of mathematical equations, but also shorter and shorter time intervals 
at which experimental data could be obtained. In the 1960's Edward Lorenz of 
MIT used a computer to numerically integrate twelve differential equations 
which described a crude model of the weather. In his simulations he found 
dramatic differences in output for subtle changes in input. [Ref. 3: p. 51] Since 
then, many other researchers have used computers to investigate chaotic 
behavior in virtually every field including mathematics, biology, chemistry, 
ecology, economics, physics, medicine, meteorology, and engineering. James 
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Gleick presents an interesting and easy-to-read development of the science of 
chaos in Chaos : Making a New Science [Ref. 4]. 

The chaotic behavior of flexible space structures has been studied by 
Moon and Li. In their analysis of a 3.5 meter long 3-D truss, they found that the 
free-play in the joints of the structure resulted in chaotic motion and modal 
frequencies considerably lower than those determined . by linear methods. 
Chaotic dynamics in this structure and other space structures might make it 
difficult to design active controls to damp out transient dynamics. [Ref. 5] In 
addition, chaotic dynamics leads to inaccuracies in linear finite element models 
which are currently used to obtain modal analysis during the design of 
spacecraft. The results of the Low Power Atmospheric Compensation 
Experiment (LACE) may lead to improvements of these finite element models. 
[Ref. 6] 

In another study by Moon and Shaw, chaotic vibration were observed in 
an elastic beam with non-linear boundary conditions. For certain values of 
forcing frequency and amplitude, the periodic motion became unstable and 
nonperiodic. The results were typical of a class of mechanical oscillators with 
play or amplitude constraining stops. [Ref. 7] 

B. GEOMETRIC AND TOPOLOGICAL METHODS OF CHAOS 

The purpose of this section is to describe several methods which have been 
developed for detecting chaotic behavior in dynamical systems. The methods 
which are discussed are time series, phase plane, frequency domain, Poincare 
section, Lyapunov exponent, and fractal correlation dimension. The first four 
methods are more qualitative in nature, and answer the question: Is chaos 
present? The last two methods, Lyapunov exponent and fractal correlation 
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dimension, are more quantitative in nature. Generally, when determining if 
chaos is present in a system, several of the above methods are used since 
applying one method by itself is not always indicative of chaotic motion. 

1. Time Domain Analysis 

The traditional way of observing signals is to view them in the time 
domain. This is simply a record of what is happening to a parameter of the 
system versus time. Sometimes the first indication of chaotic motion is that the 
parameter observed exhibits no pattern or periodicity. This test is not always 
conclusive, since the motion observed could have a long period behavior that is 
not readily observable. In addition, some nonlinear systems may appear non 
periodic, but can be broken down into several periodic signals. [Ref. 2: p. 42] 

Figure 2.1 illustrates both the nonperiodic nature of chaotic signals and 
also the sensitive dependence on initial conditions. Thompson and Stewart 
numerically integrated the Duffing's equation 

x + 0.05x + x 3 = 7.5 cos t (2.1) 

to produce the time history in the figure. The ragged appearance persists for as 
long a time as the integrations are carried out. There is never an exact 
repetition of any part of the signal— it is truly nonperiodic. The figure also 
illustrates that for slight changes in initial conditions, the trajectories quickly 
diverge. [Ref. 8: p. 4] 
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2. Frequency Domain Analysis 

Because of the problems encountered when trying to detect chaos from 
the time domain, the frequency domain, or Fourier Spectrum, is primarily used. 
Observing the signal in the frequency domain will instantly reveal if the signal is 
made up of a sum of periodic signals. Chaotic vibrations will exhibit a broad 
spectrum of frequencies in the output when the input is a single frequency 
harmonic motion. Subharmonics and Superharmonics of this frequency will 
appear not as sharp spikes like the excitation frequency but as rounded hills. 
Once again, this test is not always conclusive since multiharmonic outputs may 
result from hidden degrees of freedom of which the observer is unaware. In 
order to detect chaotic vibrations in large degree-of-freedom systems, it is 
necessary to observe changes in the spectrum while one parameter such as 
driving amplitude or frequency is varied [Ref. 2: p. 47]. 

Figures 2.2 and 2.3 show the Fourier spectrum for a periodic signal 
generated by the equation 



x = 3 cos (4k t) (2.2) 

and the chaotic signal of Equation (2.1) respectively. Notice that Figure 2.2 has 
a sharp spike only at the excitation frequency while Figure 2.3 has multiple, 

v 

broadband spikes at frequencies other than the excitation frequency of 0.16 Hz. 
3. Phase Plane 

A dynamical system consists of two parts: the state (essential 

information about a system) and the dynamic (rule that describes how the state 
evolves with time). The motion of the system can be visualized in a state space 



or phase plane. The phase plane is a construct whose coordinates are 
components of the state. The coordinates, therefore, depend on the type of 
system being studied. A mechanical system will usually use position versus 
velocity. A simple example to describe this is the motion of a pendulum. All 
that is needed to determine its motion are position and velocity. The state is 
thus a point on the phase plane. Newton's Laws provide the dynamic, 
expressed mathematically as a differential equation, that describes how the 
state evolves. Without the presence of friction, the pendulum traces out a circle 
or orbit in the phase plane. When friction is present, the orbit spirals in to a 
single point or attractor. [Ref. 3: p. 49] This describes the motion of a simple 
non chaotic system. The orbits of chaotic systems never close or retrace 
themselves. Thus, the trajectory of the orbits will tend to fill up a section of the 
phase plane. Once again, this is not always a conclusive test for chaos. 

The pseudo-phase plane is another technique which has been used 
when only one variable is able to be measured. This technique was discovered 
by scientists at the University of Santa Cruz while analyzing a dripping faucet. 
They constructed a pseudo-phase plane using a measured variable and a time- 
shifted version of the same measured variable. Takens later proved that a 
pseudo-phase plane exhibited the same characteristics as a phase plane. [Ref. 
4: p. 266] The choice of the time shift or embedding time will shrink or 
expand the orbit about the line y = x. Dvorak and Klaschka have proven that 
an embedding time between five and eleven time intervals produces a pseudo- 
phase plane with the least amount of bias [Ref. 9]. Chaotic motions exhibits the 
same characteristics in the pseudo-phase plane as they do in the phase plane. 
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The pseudo-phase plane for the periodic signal appears as a circle as 
expected (Figure 2.4). This signal constantly retraces itself. As seen in Figure 
2.5, the pseudo-phase plane for the chaotic signal of Equation (2.1) never 
retraces itself and fills up a section of the plane. 

4. Poincare Section 

The Poincare section is constructed by passing a two-dimensional 
surface through a three-dimensional phase space and observing where the 
points of the trajectory intersect the plane. This slice can be made through an 
infinite number of locations to reveal the internal structure of the trajectory. 
Periodic signals, such as an undamped pendulum, have one point on the 
Poincare section while signals with multiple harmonics produce a finite number 
of points on the Poincare section. If the Poincare section does not consist of a 
finite set of points or a closed orbit, the motion may be chaotic. In undamped or 
slightly damped chaotic systems, the Poincare section appears as a cloud of 
unorganized points. In damped chaotic systems, the Poincare section appears 
as an infinite set of highly organized points which look like parallel lines. If after 
several enlargements the structure continues to exist, the system is described 
as having a strange attractor. This embedding of structure within structure is 
referred to as a Cantor set after the mathematician George Cantor (1845 - 
1918). [Ref. 2: p. 52] 

Figure 2.6 is the Poincare section of the periodic signal. Notice that 
there is only one point in the top right hand comer. The damped, chaotic signal, 
however, appears as a set of points forming parrallel lines as shown in Figure 
2.7. If this figure were enlarged, the same pattern would emerge after each 
enlargement. This is referred to as the self-similar property. 
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5. Lyapunov Exponent 

The Lyapunov exponent is a measure of the average rate of 
divergence of trajectories initially separated by an infinitesimal amount. A 
system is chaotic if it has one or more positive Lyapunov exponents, which 
means nearby trajectories are diverging [Ref. 10: p. 323]. Lyapunov exponents 
have proven to be a useful dynamical test for chaotic systems. Since nearby 
trajectories correspond to nearly identical states, exponential divergence 
between the trajectories means that predictability of the system is rapidly lost. 
Systems with subtle, unresolvable initial differences quickly behave differently. 
[Ref. 1 1 : p. 285] The Lyapunov exponent, X, is defined as 

d(t) = d 2 Xt 
or 

^ = 7 log 2 ( 77 -) 

where : d is the initial distance between the two trajectories 

d(t) is the distance at a later time, t 

As indicated in Equation (2.4) a positive exponent means that 'd(t)' is greater 
than ’d’ and the trajectories are diverging. A zero exponent indicates 'd(t)’ and 
'd' are the same and the trajectories coincide. A negative exponent indicates 
’d(t)' is less than ’d’ and the trajectories are converging. 

A system with n degrees of freedom has a set of n exponents, called 
the Lyapunov spectrum. This spectrum characterizes the divergence of the 
trajectory in n directions. In order to understand the meaning of the Lyapunov 
spectrum, imagine a sphere of neighboring points in the phase space about a 



(2.3) 

(2.4) 
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central point on the trajectory. As the points evolve with time, this sphere will 
deform. If one axis is diverging exponentially the sphere stretches along that 
direction. If another axis is converging, the sphere shrinks along that direction. 
Algorithms are available to calculate the exponents of the spectrum from the 
measurement of a single observed state by constructing the pseudo-phase 
space. [Ref 12: p. 89] The program "Chaos" uses a method developed by 
Wolfet al [Ref. 11]. 

6. Fractal Correlation Dimension 

The fractal correlation dimension is defined as 



C(r) = r d 



(2.5) 



where: C(r) is the probability that a point on the attractor is within a 

circle, sphere, or hypersphere of radius ’r' 
r is the radius of the circle, sphere, or hypersphere 
d is the fractal correlation dimension. 

Solving for 'd' in Equation (2.5) yields 



d = lim { 



ln[C(r)] 
ln(r) J 



(2.6) 



C(r) may be calculated by constructing a sphere at each point on the 
trajectory and counting the number of points in each sphere, such that 

C(r) = lim X X H ( r - | Xi - xl ) (2.7) 

N 2 i j 
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where: 



H(s) is the Heaviside function and is equal to one if 's' is greater 
than zero and equal to zero if 's' is less than zero 
I xi - xj I is the Euclidean distance between points 
N is the total number of points [Ref. 2: p.217] 

The procedure used in the "Chaos" program by Sarigul-Klijn to 
compute the fractal correlation dimension is described below. 

1. Start with a point on the attractor and calculate the number of points 
inside a circle of radius, r. 

2. Calculate the probability C(r) by dividing this number of points by the 
total number of points in the attractor. 

3. Repeat this for several points along the attractor. 

4. Compute C(r) for several values of r. 

5. The slope of log[C(r)] versus log(r) yields d, the fractal correlation 
dimension, for the chosen number of pseudo-phase space dimensions. 
[Ref. 1: p. 88]. 

6. Repeat steps 1 through 5 for higher pseudo-phase space dimensions to 
obtain an asymptotic estimate of the fractal correlation dimension for the 
dynamical system. 

Nonchaotic systems have integer fractal correlation dimensions such as 
an equilibrium point (zero), limit cycle (one), or toroidal attractor (two). The 
fractal correlation dimension of a chaotic system is non-integer and is 
independent of the dimension of the phase space used for the calculation. [Ref. 
13: p. 101] 
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Figure 2.1 Divergence from adjacent starts. 
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III. THEORETICAL ANALYSIS 



A. DUFFING’S EQUATION 

The mathematical equation which describes the motion of a nonlinear 
elastic beam is represented by the well-known Duffing's equation. 

mv + cv + kjv + k 3 V 3 = Acos cot (3.1) 

where: m = mass of the beam 

c = damping coefficient 
ki = linear stiffness coefficient 
k 3 = nonlinear stiffness coefficient 
A = amplitude of the excitation 
(0 = frequency of the excitation 
v = transverse displacement of the beam 

Once the coefficients of the Duffing’s equation are known it can be integrated 
numerically to solve for the displacement, velocity and acceleration of the tip of 
a cantilever beam. A fourth order Runge-Kutta scheme (see Appendix A) was 
developed to perform the integration. The following methods are used to 
determine the coefficients using the properties of a beam made of 7075-76 
sheet aluminum listed in Table 3.1. 
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B. ASSUMED MODES METHOD 



In order to determine the coefficients to be used in Duffings equation, the 
assumed modes method is used. This method is an extension of the principle 
of virtual displacements and is used to produce a generalized parameter model 
that approximates the flexible behavior of the beam. A virtual displacement is 
an infinitesimal, imaginary change in the displacement of a system as shown in 
Figure 3.1. 

The dotted curve in the figure shows the virtual displacement, 5v(x,t), of the 
beam. It must satisfy the same boundary conditions as the actual displacement, 
v(x,t)--the solid curve in Figure 3.1. The virtual displacement is not a function 
of time in the same sense as the actual displacement but is rather a small 
change of configuration relative to the configuration of the beam at the time "t" 
[Ref. 14: p. 33]. 

To create the single degree of freedom model, the deflection of the beam 
may be approximated by 



v(x,t) = y(x) i}(t) (3.2) 

Similarly, the virtual displacement is given by 

8v(x,t) = 5tf(t) (3.3) 

The shape function, \|/(x), may be any admissible function which adequately 
describes the behavior of the system and possesses derivatives of an order at 
least equal to those specified by the boundary conditions. The generalized 
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displacement, d(t), is the solution to Duffing's equation. [Ref. 14: p. 34] The 
shape function which is chosen is the characteristic function for a clamped-free 
beam. 



V(4) = VI sin ‘ 4) + ex P ('f 5) + ex P ( - j (1 - £)) (3.4) 



where ^ is the beam coordinate, nondimensionalized with respect to the length, 
L, and varying between 0 and 1. [Ref 15: p. 439] 

The boundary conditions for a cantilever beam are 



BC1: 


v(0,t) = 0 


BC2: 


v'(0,t) = 0 


BC3: 


v"(L,t) = 0 


BC4: 


v'"(L,t) = 0 



(3.5) 



where (') indicates the derivative with respect to the spatial coordinate x. 

The first boundary condition states that there is no displacement at the fixed 
end while the second states that its slope must also vanish here. The third and 
fourth boundary conditions state that there must be no moment and no shear 
forces at the free end. Applying the boundary conditions of Equation (3.5) to 
the shape function of Equation (3.4) gives only a 0.7% error in amplitude at the 
midpoint. [Ref 16: p. 1014] 

The first and second derivatives of the shape function will be necessary in 
determining the coefficients for Duffing's equation. Subsequently, 



V'(5) = V2 | cos - 2.) - 1 exp ( - S- ^) + 1 exp ( - 1 (i - £)) ( 3 . 6 ) 

V"(5) = - VI sin (&£ - exp ( - 2- £) + al exp ( . TL (i . £)) (3 7 ) 
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The principle of virtual displacements states that the combined virtual work 
of all forces must be zero. More explicitly this is stated as 

8W nc - 5V +SW dmp + SW 

inert ® (3*8) 



where SW nc is the virtual work of nonconservative forces,' 8V is the change in 
strain energy of the system, 8Wdamp is the virtual work of damping forces, 
and 5Winert is the virtual work of inertia forces [Ref. 14: p. 34]. 

1. Mass Coefficient Determination 

Solving for the virtual work of inertia forces will yield the mass term, m, 
of Duffing’s Equation (3.1). Since the inertia force per unit length is distributed 
along the beam, the virtual work of inertia forces is 



8W 



inert 



■f 



pA v 8v dx 

= -pAL I xj/ 2 d£, 80 

Jo 



(3.9) 



Using the values for density, area, and length of the beam contained in Table 
3.1 and numerically integrating using MATHCAD (see Appendix B), the work 
done by inertia forces is 



SWinert = 1.0360 x 10 1 0 8t3 (3.10) 

Therefore, the mass coefficient, m, of Duffing’s equation is 0.10360 slugs. 
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2. Stiffness Coefficients Determination 

If a beam is modeled as shown in Figure 3.2 with coordinate axes x, y, 
and z and corresponding displacements u, v, and w, it can be shown [Ref. 17: p. 
22] that the strain in the x direction is given by the following equation: 



3u i 
£xx ~dx + 2 



0U 2 2 2 

^ +C-) +(^) 



dx 5x ; 



dx 



(3.11) 



Assuming that all higher order terms are negligible except for (dv/dx) 2 (small 
strain, large rotation assumption), the above equation becomes: 



e 



xA ) 2 

dx 2 dx 



(3.12) 



Applying the classical Euler-Bernoulli assumption that cross sections are 
assumed to remain planar and normal to the centroidal axis after deformation 
[Ref. 18: p. 152], it can be seen from Figure 3.3 that 



u = -y tan 0 Z 

tan 6 Z = ^ 
Ax 



(3.13) 

(3.14) 



Subsequently 



u 

3u 

dx 



- v dv 
y dx 

:. V £v 

dx 2 



(3.15) 

(3.16) 
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Substituting Equation (3.16) into Equation (3.12) yields 



f - -v 
e " y dx* 



+ U&.) 

2 dx 2 



(3.17) 



The expression for the strain energy of a beam undergoing transverse 
deflection is given by [Ref. 19: p. 169] 



V 




£xx 



dx dA 



(3.18) 



Because the beam under study is not deflected so the material passes its yield 
point, it behaves as a Hookean material such that 



Oxx 



— E e xx 



(3.19) 



Therefore the strain energy becomes 





e£ x dx dA 



(3.20) 



Squaring Equation (3.17) and substituting it into Equation (3.20) yields: 





y2(C^V) 2 

y dx 2 



V (dv_\ 2 d^v_ 
n dx ; dx 2 



4 dx 



dx dA 



(3.21) 
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If the intersection of the y and z axes lies at the centroid then 




dA = 1 



ZZ 



dA = 0 



Therefore, Equation (3.21) becomes: 



-¥/ 



V = El I (£&- ) 2 dx + EA 



dx 2 



f 



( dV ) 4 

dx 



dx 



Applying variational calculus, Equation (3.24) becomes 



5V = EI 



f d 2 vd 2 8v 

J 0 dx 2 dx 2 




/dv \ 3 dSv 
W dx 



dx 



which can be written 



f' 



5V = E I v" (8v)" dx + EA | ( V ') 3 (5v)' dx 



JO 



Substituting the derivatives of Equation (3.2) and (3.3) and 
Equations (3.6) and (3.7) into Equation (3.26) yields 

8V = El I (v") 2 $ 8£ dx + ^ I (\p') 4 i} 3 8i} dx 



(3.22) 

(3.23) 

(3.24) 



(3.25) 



(3.26) 
substituting 

(3.27) 
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Substituting in the values from Table 3.1 and numerically integrating using 
MATHCAD (See Appendix B), Equation (3.27) becomes 



8V = 1.5433 x 10- 2 8fl + 4.5943 x 10 1 fl 3 8i3 



(3.28) 



From Equation (3.28), the linear and nonlinear stiffness coefficients, kjand k 3 , 
are seen to be 0.015433 lbs/in and 45.943 lbs/in^ respectively. 

3. Damping Coefficient Determination 

If the beam is assumed to be modelled by a linear single degree of 
freedom system the damping coefficient is given by 



The values for kj and m were determined previously. In order to 
determine the viscous damping factor, the logarithmic decrement method was 
used [Ref. 14: p. 62]. The logarithmic decrement, 8, is defined by 



Cd = 2 V k i m C, 



(3.29) 



where: 



C, = viscous damping factor 




(3.30) 



where: Up = amplitude of motion at beginning of cycle 



Uq = amplitude of motion after N cycles 
N = number of cycles 
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The above equation is used to solve for the viscous damping factor by 
displacing the tip of the beam and then measuring the two amplitudes, Up and 
Uq, and counting the number of cycles between the two measurements. These 
measurements were performed on an air table in the Flexible Spacecraft 
Simulator Laboratory of NPS as shown in Figures 3.4 and 3.5. The beam was 
horizontally mounted on the table and an air pad was attached to the tip of the 
beam. The tip was then deflected and the resulting motion achieved essentially 
without friction. The values for Up, Uq, and N appear in Table 3.2. Using 
Equation (3.30), the average viscous damping factor, was determined to be 
6.67x1 O' 3 or approximately 0.7%. 

The virtual work of the damping forces is given by 



Substituting in the value for Cd obtained from Equation (3.29) and numerically 
integrating once again yields 



Therefore, the damping coefficient, c, of Equation (3.1) is 0.0025654 slugs/s. 




(3.31) 



5W damp = 2.5654 xlO'H Stf 



(3.32) 
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C. THEORETICAL RESULTS 

Appendix A contains the FORTRAN computer code used to numerically 
integrate the Duffing's equation which mathematically describes the motion of 
the flexible beam being studied. Table 3.3 summarizes the results of the 
previous section in which the coefficients of Equation (3.1) were determined. 
Equation (3.1) was then integrated using an amplitude of 3 lbs. (the force 
applied by the shaker) and a frequency of 5.6 Hz. The initial displacement and 
velocity chosen were 1 inch and 0 in/s, respectively. The results were then 
formatted to be read by the computer program "Chaos" developed by M. 
Sarigul Klijn [Ref. lj. 

Qualitative tests for chaos are satisfied by the time series, Fourier spectrum, 
pseudo-phase plane and Poincare section. Figure 3.6 shows the time history of 
the acceleration of Equation (3.1). It shows no noticeable periodic behavior and 
this is confirmed in the Fourier spectrum of Figure 3.7. The power spectrum is 
made up of many frequencies with broad band characteristics. The pseudo- 
phase plane of Figure 3.8 was obtained using an embedding time of 10. This 
embedding time, consistent with the findings of Dvorak and Klaschka [Ref 9], 
yielded the least amount of bias. Notice that the trajectory frequently crosses 
itself and fills up a section of the phase space. The Poincare section, as shown 
in Figure 3.9, is also characteristic of damped chaotic systems with its cloud of 
unorganized points filling a section of the space. It was obtained by passing a 
plane through the three-dimensional pseudo-phase space at z = 0. 

The quantitative tests for chaos were also performed. Figure 3.10 is a plot 
of the average Lyapunov exponent. This figure is obtained by time averaging 
the many computations of Equation (2.4) and then plotting them versus the 
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number of samples, 4096. The average Lyapunov exponent for this system is 
0.949. This exponent is indicative of chaotic systems. 

Figure 3.11 shows the plot from which the fractal correlation dimension is 
determined. As mentioned previously, the fractal correlation dimension is the 
asymptotic value obtained by plotting the slope of log[C(r)] versus log(r) for 
successively higher pseudo-phase space dimensions. It was determined to be 
approximately 2.4. Moon reports a value of 2.5 for his buckled beam 
experiment. [Ref. 2, p. 232]. The non-integer fractal correlation dimension is 
yet another indication of chaos in dynamic systems. 
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TABLE 3.1 Beam Properties 



Density, p 


0.101 lbm/in^ 


Modulus, E 


10.4 x lO^psi 


Length, L 


90.0 in 


Area, A 


0.24 in 2 


Thickness, t 


0.080 in 


Mom Inert, I 


0.000128 in 4 


Width, w 


3.0 in 


> 





TABLE 3.2 Logarithmic Decrement Values 



N 


Up (in) 


Uq (in) 


c 


20 


2.5625 


1.0000 


0.00749 


20 


2.3125 


0.9375 


0.00718 


20 


2.2500 


0.9375 


0.00697 


20 


2.6250 


1.0625 


0.00720 


20 


2.7500 


1.1250 


0.00711 


25 


2.2500 


0.8125 


0.00648 


25 


2.6250 


0.8750 


0.00699 


25 


2.5000 


0.8750 


0.00668 


25 


2.2500 


0.8125 


0.00648 


25 


2.2500 


0.8125 


0.00648 


30 


1.8125 


0.6250 


0.00565 


30 


2.6250 


0.7500 


0.00665 


30 


2.1875 


0.6875 


0.00614 


30 


1.9375 


0.5625 


0.00656 


30 


1.9375 


0.6250 


0.00600 
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TABLE 3.3 Coefficients of Duffing’s Equation 

m= 1.0360 x IQ' 1 
c = 2.5654 x 10-3 
kl = 1.5433 x 10- 2 
k3 = 4.5943 x 10* 




Figure 3.1 Illustration of virtual displacements 




Figure 3.2 Coordinates of a simple beam 
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y *1 x 7 *3 




Figure 3.3 Kinematics of a simple beam 
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Figure 3.4 View 1 of beam on air table 
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Figure 3.5 View 2 of beam on air table 
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time series 




POWER SPECTRUM 




Figure 3.6 Time series of Duffing's 
Equation. 



Figure 3.7 Fourier spectrum 
Duffing's equation. 




Figure 3.8 Pseudo-phase plane of Figure 3.9 Poincare section 
Duffing's Equation. of Duffing's equation. 
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LYAPUNOV EXPONENT 




Figure 3.10 Lyapunov Exponent of Duffing's equation. 




Figure 3.11 



Fractal Correlation Dimension of Duffing's equation 
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IV. EXPERIMENTAL ANALYSIS 



A. EXPERIMENTAL SET-UP 

The aluminum beam with properties listed in Table 3.1 is instrumented with 
strain gauges and an accelerometer as shown in Figure 4.1. The type of strain 
gauges used are one-quarter inch, 350 ohm, straight geometry gauges from the 
Micro Measurements Group. The accelerometer is a Sensotec model JTF flat 
pack. This accelerometer has a weight of one ounce, a peak range of ±5 g's, a 
sensitivity of 10 mV/g, and a usable frequency range of DC to 150 Hz. 

The beam is mounted vertically to an MB Dynamics Model PM500A 
Vibration Exciter. The outputs from the strain gauges and accelerometer are 
then fed into a bridge balancing box and amplifier. The bridge box is capable of 
handling four inputs, one from the accelerometer and three from the strain 
gauges. From the bridge box, the signal goes to an analog-to-digital converter. 
From there the acceleration and strain data is stored on a Zenith computer. 
Mechanical stops placed 58 inches up from the base prevent the beam from 
buckling past one and three-fourths inches on either side. (Figure 4.2) 

B. EXPERIMENTAL PROCEDURE 

The accelerometer is placed at location A1 in Figure 4.1. The 
accelerometer and the three strain gauges are then connected to the data 
acquisition system. The beam is held in the upright vertical position in order to 
zero out all the readings from the gauges and the accelerometer using the 
bridge balancing circuit. The beam is then released while simultaneously 
starting the exciter. The exciter provides an excitation frequency of 5.6 Hz and 
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a displacement amplitude of one-quarter inch peak-to-peak. This procedure is 
then repeated twice, placing the accelerometer at locations A2 and A3. Strain 
gauge measurements at location SI, S2, and S3 are obtained for each run. 

C. EXPERIMENTAL RESULTS 

The data from each of the runs is transferred into a format readable by the 
"Chaos" program. It is then analyzed using the techniques discussed in 
Chapter II. These results appear at the end of this chapter. As shown in 
Figures 4.3 through 4.10, the four qualitative tests for chaos seem to be met: 

1. Erratic and unpredictable time histories from both accelerometers and 
strain gauges. 

2. Multiple broadbanded power spectrum spikes. 

3. Crossing trajectories and filled up pseudo-phase space. 

4. A cloud of unorganized points in the Poincare section characteristic of 
lightly damped systems. 

These results appear similar to the chaotic results obtained by Sarigul-Klijn 
[Ref. 1] in his analysis of helicopter flight test data using higher harmonic 
control. 

The quantitative tests show that the location of the type of sensor chosen 
impacts whether or not chaotic motion is found. Strain gauges only indicate 
chaos if placed near the root, while accelerometers only indicate chaos if 
placed near the tip. Figures 4.11 through 4.14 are plots of the average 
Lyapunov exponents and fractal correlation dimensions. A summary of the 
Lyapunov exponents and fractal correlation dimensions for each sensor appear 
in Tables 4.1 and 4.2. The Lyapunov exponents are all positive for the strain 
gauges. For the accelerometers, the Lyapunov exponent diminishes as the 
accelerometer is placed closer to the root. This indicates that measuring 
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accelerations near the root reduces the chances of obtaining chaotic results. 
This is anticipated since the beam does not experience large deflections in this 
area and can be characterized by linear equations. Likewise, the fractal 
correlation dimensions of the sensors indicate that accelerometers should be 
placed closer to the tip and strain gauges closer to the root. The accelerometer 
nearest the root has an integer fractal correlation dimension, like that of a non- 
chaotic system. The two nearer to the tip have non-integer fractal dimensions, 
characteristic of chaotic systems. Conversely, the two strain gauges near the 
root have non-integer fractal dimensions, while the one closer to the tip has an 
integer fractal dimension. 
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TABLE 4.1 Lyapunov exponents 



SI 


0.272 


A1 


0.364 


S2 


0.237 


A2 


0.0397 


S3 


0.390 


A3 


-0.0282 



TABLE 4.2 Fractal Correlation Dimensions 



SI 


4.8 


A1 


4.4 


S2 


4.4 


A2 


4.4 


S3 


4.0 


A3 


4.0 
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TIP 



ROOT 



( 90" 1 



□ □ □ □ □ □ 



A1 S3 A2 A3 S2 SI 



Strain Gauges 

\* a i 

| — 14.5" — | 

I 66" 1 

Accelerometer 

[2.Zi " 

| 30.5" 1 

I 5 4" H 

Figure 4.1 Sensor Locations 




Figure 4.2 Experimental Set-up 
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(C) 

Figure 4.3 Time series of strain gauges 

(a) SI (b) S2 (c) S3 
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TIME SERIES TIME SERIES 




(a) 



(b) 



TIME SERIES 




(c) 

Figure 4.4 Time series of accelerometers 

(a) A1 (b) A2 (c) A3 
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Figure 4.5 Pseudo-phase planes of strain gauges 

(a) SI (b) S2 (c) S3 
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Phase plane 



phase Plane 




e 



(b) 

phase plane 



-b - 

-e - 
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PiTl 

(C) 

Figure 4.6 Pseudo-phase planes of accelerometers 

(a) A1 (b) A2 (c) A3 
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POWER 



POWER SPECTRUM POWER SPECTRUM 




(a) 



(b) 



POWER SPECTRUM 




(c) 

Figure 4.7 Fourier spectrums of strain gauges 

(a) SI (b) S2 (c) S3 
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POufP 



POWER SPECTRUM 



POWER SPECTRUM 




(a) 



(b) 



POWER SPECTRUM 




(c) 



Figure 4.8 Fourier spectrums of accelerometers 

(a) A1 (b) A2 (c) A3 
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POINCARE SECTION Of 3-D SPACE 




POINCARE SECTION OF 3-D SPACE 




(b) 



POINCARE SECTION OF 3-D SPACE 




(c) 

Figure 4.9 Poincare sections of strain gauges 

(a) SI (b) S2 (c) S3 
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F i f - TAUI 



POINCARE SECTION OP 3-0 SPACE 



POINCARE SECTION OF 3-D SPACE 
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(b) 



POINCARE SECTION OF 3-D SPACE 




(c) 

Figure 4.10 Poincare sections of accelerometers 

(a) A1 (b) A2 (c) A3 
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trAPUNOV EXPONENT 



lyapunov exponent 





(a) 



(b) 



LYAPUNOV EXPONENT 




(c) 

Figure 4.11 Average Lyapunov exponents of strain gauges 

(a) SI (b) S2 (c) S3 
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ltapunov exponent 




LYAPUNOV EXPONENT 




(c) 

Figure 4.12 Average Lyapunov exponents of accelerometers 

(a) A1 (b) A2 (c) A3 



51 





(a) (b) 




Figure 4.13 



(c) 

Fractal correlation dimension of strain gauges 

(a) SI (b) S2 (c) S3 
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a d h 





(a) (b) 




(c) 

Figure 4.14 Fractal correlation dimension of accelerometers 

(a) A1 (b) A2 (c) A3 
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V. CONCLUSIONS AND SCOPE FOR FUTURE 

RESEARCH 



For the first time geometric and topological methods of chaos were applied 
to data obtained simultaneously by accelerometers and strain gauges. It was 
determined that data from both accelerometers and strain gauges provided 
similar indications to the presence of chaos when analyzed using time series, 
Fourier spectrums, pseudo-phase planes, and Poincare sections. Using these 
techniques, both types of sensors indicated chaos without regard to the location 
of the sensor. However, the Lyapunov exponents and fractal correlation 
dimensions indicated that chaotic motion is only detected by accelerometers if 
they are placed near the tip and that chaotic motion is only detected by strain 
gauges if the are placed near the root. 

The presence of chaos in highly flexible spacecraft structures was also 
investigated. It was determined that the geometric nonlinearities of these 
structures may result in chaotic vibrations and dynamic instabilities. These 
instabilities may then affect control mechanisms of the spacecraft. Further 
analysis could be performed in using the geometric and topological methods of 
chaos to analyze the effectiveness of these mechanisms to adequately control 
the spacecraft and provide improvements into the way the control is performed. 
Chaotic methods may also have an application in updating the linear finite 
elements models of flexible space structures with more accurate nonlinear 
models. 
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APPENDIX A 



TITLE t OUirXMO 

DATE i 07 DEC 1001 



AUTHOR; ROBERT Q. VAUGHAN 

SYSTEM! VAX 

COMPILER I HASP 



DEE CRI FT 1 091 1 

This program solves Duffing'* equation 



s*d~2/dt"2(x) ♦ c*d/dt(x) ♦ ll*x ♦ k3*x~3 » Acos(v*t) 



by transforming it into two ordinary differential Equation* 
d/dt(x) - y 



d/dt(y) m acoe( v* t ) - c*y - kl*x - *3*x“3 

and using a fourth order Stunga— Rutta Eoktac to intagrata to 
find tha displacement, x, and th* velocity, y. 

Thi amplitude and frequency of th# forcing function, tha atap 
Elsa for int*gration r the number of lntagratlon step*, and tha 
initial valuas for displacement and velocity ara antarad in 
that ordar in a data file titled "DUFFIfeflT. 1 The data file 
titled "Durrour” contain* the integrated valuat for displacamant 
and velocity In a format compatible with the program "CHAOS" 
by n» florigul-Rli jn. 

To alter thi* program to aolva an aquation other than Duffing'* 
aboTc, the parameter, NS, must be changed to correspond to the 
number of 0DE r e and the subroutine •DIFFIQ" must be changed 
to Include the proper equations. 



Two subroutine* are employed. 

DIFFIQ calculates tha valuta for tha ODE'* 

RK4 is the Runga-Kutta routine 

A* written the program do#* not usa an adaptive etep~aise 

The variable a ara defined aa follows t 

Ni - the number of ODE' a in the system 

NS * the number of steps (iterations) to be performed 

B - the desired stepeise 

A - amplitude of the forcing function 

F - frequency of tha forcing function 

TIME, T - time 

I, INDEX - integer counters 

X - an array corresponding to the Independent variables, 
i.e. X(l)-X, X(2)-T, etc. 
root - an array cor responding to tha OW'f 



PROGRAM DUFFING 
PARAMETER (NE-2) 

IMPLICIT DOUBLE PBICI£IOM( A-H,X,M, 0~Z ) 
CHARACTER *7 DUMMY 
DIMENSION X(NE) ,XD0T(NE) 
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0PEN(7, riLE-'DUrriNIT' , STATUS- ' OLD ' ) 
OPEN! 6 , FILE— ’ DUTPOUT* , STATUS- ' OLD ' ) 



TIME-0. 0D0 
DUMMY-' A1234A' 

READ ( 7 , * ) A,F,B,NS 
DO 5 I-1,HE 

READ( 7 , * ) X(I) 

5 CONTINUE 

WRITE! 6,6) A,T 
WRITB(8,7) B,NE 
WRITE! 8,8) X(1),X(2) 

6 FORMAT (IX, 'FORCING FUNCT AMPLITUDE - ' , F4 . 1 , 5X, ' FREQ - ',F6.3) 

7 FORMAT (IX, 'STEF SIZE - ' , F6 . 4 , 5X, ' NUMBER Or 8TBP8 - ',15) 

6 FORMAT ( IX, ' INXT DI6P - ' , P4 . 1 , SX, ' IKIT VEL - ',F4.1> 

CALL DIFFEQ(X,XDOT,A,F,TIME) 

WRITE (8, 100) DUMMY , TIME,X( 1 ) , X( 2 ) , XDOT( 2 ) 

DO 15 I— 1 , NS 

DO 10 INDEX-1,4 

CALL DIFFEQ(X,XDOT,A,F,TIKE) 

CALL RK4(TIME,X,XDOT,H, INDEX) 

10 CONTINUE 

WRITE ( 8 , 100 ) DUMMY, TIME, X( 1 ) ,X( 2 ) ,XDOT( 2 ) 

15 CONTINUE 

100 FORMAT (A7, F6.3, Ell. 4, 2E12.4) 

STOP 

END 

************************************************************************ 
SUBROUTINE RK4(T,X,XD0T,H, INDEX) 

* This subroutine performs numerical integration using 4th order 

• Runge-Kutta methods. It vie developed by I.M. Robb. 

PARAMETER! NE- 2) 

IMPLICIT DOUBLE PRECISION! A-H, K,M, O-Z ) 

DIMENSION X(NE) ,XDOT(NE) , SAVED ( NE ) ,SAVEX(NE) 

GO TO (1.2.3. 4) .INDEX 

1 DO 10 I-1,NE 

3AVEX< I )-X( I ) 

SAVED! I ) — XDOT ( I ) 

10 X(I)-SAVEX(I)+0.5D0*H*XDOT(I) 

T-T+0 . 5D0*H 
RETURN 
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2 DO 20 I-1,NE 

SAVED ( I ) —SAVED ( X ) + 2 . ODO*XDOT( I ) 
20 X( I )«SAVEX( I ) +0 . 5DO*H*XDOT( I ) 

RETURN 

3 DO 30 I«1,NE 

SAVED ( I ) -SAVED ( I )-r2 . ODO*XDOT< I ) 
30 X( I )-£AVEX( I ) +H*XDOT ( I ) 

T«T+0.5D0*H 

RETURN 

4 DO 40 I-1,NE 

4 0 X ( I ) -SAVEX ( I ) +B/S . 0D0 * ( SAVED ( I ) +XDOT ( I ) ) 

RETURN 

END 



SUBROUTINE DIFFEQ(X,XDOT, A, F , T) 

This subroutine obtains values for the ode's for Duffing' s 
equation and must be changed for a new system of equations. 

PARAMETER (NE-2) 

IMPLICIT DOUBLE PRECISION A-H , K, M # 0-Z ) 

DIMENSION X( NE ) , XD0T( NE ) 

PI-4 . 0D0*DATAN( 1 . 0D0 ) 

Be Careful of units! Kl - LBf/in. K3 - LBf/in*3 

M - slugs C - slugs/s 

M-l . 0360D-01 
C-2 . 56 5 4D-03 
Kl-1 . 54 33D-02 
K3-4 . 59 43D+01 

Divide M by 12 in/ft so that XD0T(2) has units in/s~2 

M-M/12.0D0 

XDOT ( 1 ) —X ( 2 ) 

XDOT(2MA*COS(F*2.0DO*PI*T)-C/12.0DO*X(2)-K1*X(1)-K3*X(1)**3)/M 

RETURN 

END 
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APPENDIX B 
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